//////////////////////////////////////////////////////////////////////
//                                                                  //
//                                                                  //
//              Copyright (C), 2000  SEIKO EPSON Corp.              //
//              ALL RIGHTS RESERVED                                 //
//                                                                  //
//                                                                  //
//   file name : log.c                                              //
//                                                                  //
//   Revision history                                               //
//       00/02/10    first release                 M.Igarashi       //
//                                                                  //
//////////////////////////////////////////////////////////////////////

#include <math.h>
#include <smcvals.h>

#define	K0	-0.18327870372215593212e2
#define	K1	 0.934639006428585382474e2
#define	K2	-0.90174691662040536328986e2
#define	S0	-0.20733487895513939345e2
#define	S1	 0.61761065598471302843e2
#define	S2	-0.45087345831020305748486e2



// ALGORITHM                
//1. Argument Reduction: find k and xfrac such that 
//			x = 2^k * (xfrac), 
//	   where  sqrt(2)/2 < xfrac < sqrt(2) .
//
//2. Approximation of log(xfrac).
//	Let x = xfrac/(1+xfrac) ; based on log(xfrac) = log(1+x) - log(1-x)
//		 = 2s + 2/3 x**3 + 2/5 x**5 + .....,
//	     	 = 2s + x*R
//	2x = xfrac-1 - x*(xfrac-1) = (xfrac-1) - hfsq + s*hfsq, where hfsq = (xfrac-1)*(xfrac-1)/2.
//	In order to guarantee error in log below 1ulp, we compute log by
//		log(xfrac) = xfrac-1 - (hfsq - s*(hfsq+R)).
//	
//3. Finally,  log(x) = k*ln2 + log(xfrac).  
//

//  the architecture of double floating point
//
//   0 1          11                  31 32                               63 bit
//   -----------------------------------------------------------------------
//  | |   exponent  |                    fraction                           |
//   -----------------------------------------------------------------------
//
//  |               |                   |                                   |
//  |    12 bits          20 bits       |             32 bits               |
//  |            lower word             |            higher word            |
//
//         bit    0         sign bit         (  1 bit  )
//              1 - 11      exponent part    ( 11 bits )
//             12 - 63      fraction part    ( 52 bits )
//
//

double log(double dfX){

	DBL_LNG       ckarg;       /* exchange the type to check the argument   */
	long lX;
	unsigned long ulHx;
	double	dfHfsq;
	double	dfXfrac,dfRet,dfTemp;
	int	iK;
	
	GETW_L(lX,dfX);		// get low
	GETW_H(ulHx,dfX);		// get high
	
/* check the argument */
	 ckarg._D = dfX;
	 CHECK_ARG( &ckarg );                     /* check NaN */
	                                          /* error: error code 33 domain error */
	if ( ckarg.st._LL == P_INF.st._LL ) {    /* x = Inf */
	    errno = ERANGE;                       /* 34: range error */
	    return( P_INF._D );
	} else if ( ( ckarg.st._LL == N_INF.st._LL ) || ( (lX&0x80000000) !=0 ) ) {  /* x < 0 */
	    errno = EDOM;                         /* 33: domain error */
	    return( NAN._D );
	} else if ( lX ==0 &&  ulHx == 0.0 ) {
	    errno = ERANGE;                       /* 34: range error */
	    return( N_INF._D );
	}
	dfX = frexp( dfX, &iK );

	GETW_L(lX,dfX);		// get low
	GETW_H(ulHx,dfX);		// get high
	
	if( (lX&0x80000000) !=0 || lX < 0x3ff6a09e || ((lX == 0x3ff6a09e) && (ulHx < 0x667f3bcd) )){
		dfX=2*dfX;
		iK=iK-1;
	}

	// calc log
	dfXfrac=(dfX-1)/(dfX+1);
	dfHfsq=dfXfrac*dfXfrac;
	dfRet=K0*dfHfsq+K1;
	dfRet=dfRet*dfHfsq+K2;
	dfTemp=dfHfsq+S0;
	dfTemp=dfTemp*dfHfsq+S1;
	dfTemp=dfTemp*dfHfsq+S2;
	dfRet=dfRet/dfTemp;
	dfRet*=dfXfrac;
	dfRet+=iK*LSS;
	return dfRet;       
}
